#include "theor_bispectrum_5.h"

double *read_transfer(char filename[], int first){
  double *ptf_arr;
  size_t n1, n2, n3;
  FILE *fp;
  int npol, nl, np, i;
  int *larr;
  double *parr;

  fp = fopen(filename, "rb");
  if (fp == NULL) {
    fprintf(stderr, "Cannot open file %s\n", filename);
    exit(1);
  }
  n1 = fread(&npol, sizeof(int), 1, fp);
  n2 = fread(&nl, sizeof(int), 1, fp);
  n3 = fread(&np, sizeof(int), 1, fp);
  if (n1 != 1 || n2 != 1 || n3 != 1) {
    fprintf(stderr, "Read error: %s\n", filename);
    fclose(fp);
    exit(1);
  }
  //printf("%d\t%d\t%d\n", npol, nl, np);
  if (first){
    tfvars.num_pol = npol;
    tfvars.num_l = nl;
    tfvars.num_p = np;
  }
  else
    if (tfvars.num_pol != npol || tfvars.num_l != nl || tfvars.num_p != np) {
      fprintf(stderr, "Error: adiabatic and isocurvature transfer functions should have same polarizations, l-values and p-values.\n");
      fclose(fp);
      exit(1);
    }
  larr = (int*) malloc(tfvars.num_l * sizeof(int));
  parr = (double*) malloc(tfvars.num_p * sizeof(double));
  ptf_arr = (double*) malloc(tfvars.num_pol*tfvars.num_l*tfvars.num_p * sizeof(double));
  if (larr == NULL || parr == NULL || ptf_arr == NULL){
    fprintf(stderr, "Memory allocation failed.\n");
    fclose(fp);
    exit(1);
  }
  n1 = fread(larr, sizeof(int), tfvars.num_l, fp);
  n2 = fread(parr, sizeof(double), tfvars.num_p, fp);
  n3 = fread(ptf_arr, sizeof(double), tfvars.num_pol*tfvars.num_l*tfvars.num_p, fp);
  fclose(fp);
  if (n1 != tfvars.num_l || n2 != tfvars.num_p || n3 != tfvars.num_pol*tfvars.num_l*tfvars.num_p) {
    fprintf(stderr, "Read error: %s\n", filename);
    exit(1);
  }
  if (first){
    tfvars.l_arr = larr;
    tfvars.p_arr = parr;
  }
  else{
    for (i = 0; i < tfvars.num_l; i++)
      if (tfvars.l_arr[i] != larr[i]){
	fprintf(stderr, "Error: adiabatic and isocurvature transfer functions should have same polarizations, l-values and p-values.\n");
	exit(1);
      }
    for (i = 0; i < tfvars.num_p; i++)
      if (tfvars.p_arr[i] != parr[i]){
	fprintf(stderr, "Error: adiabatic and isocurvature transfer functions should have same polarizations, l-values and p-values.\n");
	exit(1);
      }
    free(larr);
    free(parr);
  }
  //for (i = 0; i < tfvars.num_l; i++)
  //  printf("%d\n", larr[i]);
  //for (i = 0; i < tfvars.num_p; i++)
  //  printf("%g\n", parr[i]);
  return ptf_arr;
}

void make_splines(double *ptf_arr){
  int i, j, k, l;
  double tf_arr[tfvars.num_l][tfvars.num_pol][tfvars.num_iso][tfvars.num_p];

  tfvars.spline_arr = (gsl_spline**) malloc(tfvars.num_iso*tfvars.num_pol*tfvars.num_l * sizeof(gsl_spline*));
  tfvars.acc_arr = (gsl_interp_accel**) malloc(tfvars.num_iso*tfvars.num_pol*tfvars.num_l * sizeof(gsl_interp_accel*));
  if (tfvars.spline_arr == NULL || tfvars.acc_arr == NULL){
    fprintf(stderr, "Memory allocation failed.\n");
    exit(1);
  }
  for (k = 0; k < tfvars.num_iso; k++)
    for (j = 0; j < tfvars.num_pol; j++)
      if (tfvars.iso_summ[k] && tfvars.pol_summ[j])
	for (i = 0; i < tfvars.num_l; i++)
	  for (l = 0; l < tfvars.num_p; l++)
	    tf_arr[i][j][k][l] = *(ptf_arr + j + i*tfvars.num_pol
				   + l*tfvars.num_l*tfvars.num_pol
				   +k*tfvars.num_p*tfvars.num_l*tfvars.num_pol);
  //printf("%g\t%g\n", tf_arr[4][0][0][1000], tf_arr[4][1][0][1000]);
  free(ptf_arr);
  for (k = 0; k < tfvars.num_iso; k++)
    for (j = 0; j < tfvars.num_pol; j++)
      if (tfvars.iso_summ[k] && tfvars.pol_summ[j])
	for (i = 0; i < tfvars.num_l; i++){
	  tfvars.acc_arr[k*tfvars.num_pol*tfvars.num_l + j*tfvars.num_l + i]
	    = gsl_interp_accel_alloc();
	  tfvars.spline_arr[k*tfvars.num_pol*tfvars.num_l + j*tfvars.num_l + i]
	    = gsl_spline_alloc(gsl_interp_cspline, tfvars.num_p);
	  gsl_spline_init(tfvars.spline_arr[k*tfvars.num_pol*tfvars.num_l
					    + j*tfvars.num_l + i],
			  tfvars.p_arr, tf_arr[i][j][k], tfvars.num_p);
	}
}

void free_splines(void){
  int i, j, k;

  for (k = 0; k < tfvars.num_iso; k++)
    for (j = 0; j < tfvars.num_pol; j++)
      if (tfvars.iso_summ[k] && tfvars.pol_summ[j])
	for (i = 0; i < tfvars.num_l; i++){
	  gsl_spline_free(tfvars.spline_arr[k*tfvars.num_pol*tfvars.num_l
					    + j*tfvars.num_l + i]);
	  gsl_interp_accel_free(tfvars.acc_arr[k*tfvars.num_pol*tfvars.num_l
					       + j*tfvars.num_l + i]);
	}
  free(tfvars.spline_arr);
  free(tfvars.acc_arr);
}

void init_bess_table(double x_max, double z_spacing){
  double k_max = tfvars.p_arr[tfvars.num_p-1];
  double z_max = k_max * x_max;
  int num_z = ceil(z_max/z_spacing);
  init_spline_bessel(tfvars.num_l, tfvars.l_arr, num_z, z_max);
}

void init_theor_bispec(char **filename){
  int i, first, offset;
  double *ptf_arr, *ptf_tot_arr;

  first = 1;
  for (i = 0; i < tfvars.num_iso; i++)
    if (tfvars.iso_summ[i]){
      ptf_arr = read_transfer(filename[i], first);
      if (first){
	offset = tfvars.num_pol * tfvars.num_l * tfvars.num_p;
	ptf_tot_arr = (double*) malloc(tfvars.num_iso*offset * sizeof(double));
	if (ptf_tot_arr == NULL){
	  fprintf(stderr, "Memory allocation failed.\n");
	  exit(1);
	}
	first = 0;
      }
      memcpy(ptf_tot_arr + i*offset, ptf_arr, offset * sizeof(double));
      free(ptf_arr);
    }
  make_splines(ptf_tot_arr);
  init_bess_table(XMAX, Z_SPACING);
}

void free_all_splines(void){
  free_splines();
  free_bessel_splines();
}

double ki_integrand(double k, void *pparams){
  int pol, iso;
  double tf, bf;
  funcparams params;

  params = *(funcparams *)pparams;
  pol = params.pol[0];
  iso = params.iso[0];
  //bf = sph_bess_func(tfvars.l_arr[params.li[0]], k * params.x);
  bf = spline_bessel(params.li[0], k * params.x);
  //tf = - spline_bessel(params.li[index], k * 14000.) / 3.;
  tf = tfvars.tf_corr[iso]
    * gsl_spline_eval(tfvars.spline_arr[iso*tfvars.num_pol*tfvars.num_l
					+ pol*tfvars.num_l + params.li[0]], k,
		      tfvars.acc_arr[iso*tfvars.num_pol*tfvars.num_l +
				     pol*tfvars.num_l + params.li[0]]);
  // Note CAMB uses sign convention where SW limit is -j_l/3,
  // which means Delta T/T = -1/3 Phi.
  // Moreover, CAMB's transfer function is for R, not Phi, which means
  // that we have to multiply by 5/3 to get the one for Phi.
  // While for the isocurvature case the correction factor is 5.
  return tf * bf * pow(k/tfvars.pivot, params.power);
}

double tf_integrand(double k, void *pparams){
  int pol1, pol2, iso1, iso2;
  double tf1, tf2;
  funcparams params;

  params = *(funcparams *)pparams;
  pol1 = params.pol[0]; pol2 = params.pol[1];
  iso1 = params.iso[0]; iso2 = params.iso[1];
  tf1 = tfvars.tf_corr[iso1]
    * gsl_spline_eval(tfvars.spline_arr[iso1*tfvars.num_pol*tfvars.num_l
					+ pol1*tfvars.num_l + params.li[0]], k,
		      tfvars.acc_arr[iso1*tfvars.num_pol*tfvars.num_l
				     + pol1*tfvars.num_l + params.li[0]]);
  // CAMB's transfer function is for R, not Phi, which means
  // that we have to multiply by 5/3 to get the one for Phi.
  // While for the isocurvature case the correction factor is 5.
  if (pol1 == pol2 && iso1 == iso2)
    tf2 = tf1;
  else
    tf2 = tfvars.tf_corr[iso2]
      * gsl_spline_eval(tfvars.spline_arr[iso2*tfvars.num_pol*tfvars.num_l
					  + pol2*tfvars.num_l+params.li[0]], k,
			tfvars.acc_arr[iso2*tfvars.num_pol*tfvars.num_l
				       + pol2*tfvars.num_l + params.li[0]]);
  return tf1 * tf2 * pow(k/tfvars.pivot, params.power);
}

void init_kintvars_x(void){
  double *x1, *x2, *x3;

  kintvars.num_x = PTS_MIN*NSUB_MIN + PTS_MID*NSUB_MID + PTS_MAX*NSUB_MAX;
  kintvars.x_arr = (double*) malloc(kintvars.num_x * sizeof(double));
  if (kintvars.x_arr == NULL){
    fprintf(stderr, "Memory allocation failed.\n");
    exit(1);
  }
  x1 = get_int_points(XMIN, XMID1, PTS_MIN, NSUB_MIN);
  x2 = get_int_points(XMID1, XMID2, PTS_MID, NSUB_MID);
  x3 = get_int_points(XMID2, XMAX, PTS_MAX, NSUB_MAX);
  memcpy(kintvars.x_arr, x3, PTS_MAX*NSUB_MAX*sizeof(double));
  memcpy(kintvars.x_arr+PTS_MAX*NSUB_MAX, x2, PTS_MID*NSUB_MID*sizeof(double));
  memcpy(kintvars.x_arr+PTS_MAX*NSUB_MAX+PTS_MID*NSUB_MID, x1,
	 PTS_MIN*NSUB_MIN*sizeof(double));
  free(x1); free(x2); free(x3);
}

void init_kintvars_w(void){
  int i;
  double *Wi_min, *Wi_mid, *Wi_max;
  double width;

  Wi_min = get_weights(PTS_MIN);
  width = (XMID1 - XMIN) / NSUB_MIN;
  for (i = 0; i < PTS_MIN; i++)
    Wi_min[i] *= width/2.;
  Wi_mid = get_weights(PTS_MID);
  width = (XMID2 - XMID1) / NSUB_MID;
  for (i = 0; i < PTS_MID; i++)
    Wi_mid[i] *= width/2.;
  Wi_max = get_weights(PTS_MAX);
  width = (XMAX - XMID2) / NSUB_MAX;
  for (i = 0; i < PTS_MAX; i++)
    Wi_max[i] *= width/2.;

  kintvars.w_arr = (double*) malloc(kintvars.num_x * sizeof(double));
  if (kintvars.w_arr == NULL){
    fprintf(stderr, "Memory allocation failed.\n");
    exit(1);
  }
  for (i = 0; i < NSUB_MAX; i++)
    memcpy(kintvars.w_arr+i*PTS_MAX, Wi_max, PTS_MAX*sizeof(double));
  for (i = 0; i < NSUB_MID; i++)
    memcpy(kintvars.w_arr+PTS_MAX*NSUB_MAX+i*PTS_MID, Wi_mid,
	   PTS_MID*sizeof(double));
  for (i = 0; i < NSUB_MIN; i++)
    memcpy(kintvars.w_arr+PTS_MAX*NSUB_MAX+PTS_MID*NSUB_MID+i*PTS_MIN, Wi_min,
	   PTS_MIN*sizeof(double));
  free(Wi_min); free(Wi_mid); free(Wi_max);
}

int adaptive_step(int li, int stepstep, int stepmax){
  int step;

  step = li/stepstep + 1;
  if (step > stepmax)
    step = stepmax;
  //printf("%d\t%d\n", li, step);
  return step;
}

void init_kintvars_l(int limin, int limax, int stepstep, int stepmax,
		     int do_shapes){
  int li1, count;

  kintvars.num_li = 0;
  for (li1 = limin; li1 <= limax; li1 += adaptive_step(li1,stepstep,stepmax))
    kintvars.num_li++;
  kintvars.li_arr = (int*) malloc(kintvars.num_li * sizeof(int));
  if (kintvars.li_arr == NULL){
    fprintf(stderr, "Memory allocation failed.\n");
    exit(1);
  }
  count = 0;
  for (li1 = limin; li1 <= limax; li1 += adaptive_step(li1,stepstep,stepmax)){
    kintvars.li_arr[count] = li1;
    count++;
  }
  kintvars.alpha_arr = (double*) malloc(tfvars.num_iso * tfvars.num_pol
                       * kintvars.num_li * kintvars.num_x * sizeof(double));
  kintvars.beta_adi_arr = (double*) malloc(tfvars.num_iso * tfvars.num_pol
             		   * kintvars.num_li * kintvars.num_x * sizeof(double));
  kintvars.beta_iso_arr = (double*) malloc(tfvars.num_iso * tfvars.num_pol
             		   * kintvars.num_li * kintvars.num_x * sizeof(double));
  if (kintvars.alpha_arr == NULL || kintvars.beta_adi_arr == NULL || kintvars.beta_iso_arr == NULL){
    fprintf(stderr, "Memory allocation failed.\n");
    exit(1);
  }
  if (do_shapes > 1){
    kintvars.gamma_arr = (double*) malloc(tfvars.num_iso * tfvars.num_pol
			   * kintvars.num_li * kintvars.num_x * sizeof(double));
    kintvars.delta_arr = (double*) malloc(tfvars.num_iso * tfvars.num_pol
             		   * kintvars.num_li * kintvars.num_x * sizeof(double));
    if (kintvars.gamma_arr == NULL || kintvars.delta_arr == NULL){
      fprintf(stderr, "Memory allocation failed.\n");
      exit(1);
    }
  }
}

void calculate_k_integrals(funcparams *pparams, int do_shapes){
  int lii, xi, i, j;

  for (j = 0; j < tfvars.num_iso; j++)
    for (i = 0; i < tfvars.num_pol; i++)
      if (tfvars.iso_summ[j] && tfvars.pol_summ[i]){
	pparams->pol[0] = i;
	pparams->iso[0] = j;
	for (lii = 0; lii < kintvars.num_li; lii++){
	  pparams->li[0] = kintvars.li_arr[lii];
	  for (xi = 0; xi < kintvars.num_x; xi++){
	    pparams->x = kintvars.x_arr[xi];
	    pparams->power = 2.;
	    kintvars.alpha_arr[j*tfvars.num_pol*kintvars.num_li*kintvars.num_x
			       + i*kintvars.num_li*kintvars.num_x
			       + lii*kintvars.num_x + xi]
	      = int_1d(&ki_integrand, pparams->qmin, pparams->qmax, PTS_KINT,
		       NSUB_KINT, pparams);
	    pparams->power = -1. + tfvars.red_tilt_adi;
	    kintvars.beta_adi_arr[j*tfvars.num_pol*kintvars.num_li*kintvars.num_x
			      + i*kintvars.num_li*kintvars.num_x
			      + lii*kintvars.num_x + xi]
	      = int_1d(&ki_integrand, pparams->qmin, pparams->qmax, PTS_KINT,
		       NSUB_KINT, pparams) * tfvars.rel_ampl[j][j];
            pparams->power = -1. + tfvars.red_tilt_iso;
	    kintvars.beta_iso_arr[j*tfvars.num_pol*kintvars.num_li*kintvars.num_x
			      + i*kintvars.num_li*kintvars.num_x
			      + lii*kintvars.num_x + xi]
	      = int_1d(&ki_integrand, pparams->qmin, pparams->qmax, PTS_KINT,
		       NSUB_KINT, pparams) * tfvars.rel_ampl[j][j] * tfvars.rel_ampl_iso;
	    if (do_shapes > 1 && tfvars.iso_summ_nl[j] &&tfvars.pol_summ_nl[i]){
	      pparams->power = 1. + tfvars.red_tilt[j][j]/3.;
	      kintvars.gamma_arr[j*tfvars.num_pol*kintvars.num_li*kintvars.num_x
				 + i*kintvars.num_li*kintvars.num_x
				 + lii*kintvars.num_x + xi]
		= int_1d(&ki_integrand, pparams->qmin, pparams->qmax, PTS_KINT,
			 NSUB_KINT, pparams) * pow(tfvars.rel_ampl[j][j],1./3.);
	      pparams->power = 2.*tfvars.red_tilt[j][j]/3.;
	      kintvars.delta_arr[j*tfvars.num_pol*kintvars.num_li*kintvars.num_x
				 + i*kintvars.num_li*kintvars.num_x
				 + lii*kintvars.num_x + xi]
		= int_1d(&ki_integrand, pparams->qmin, pparams->qmax, PTS_KINT,
			 NSUB_KINT, pparams) * pow(tfvars.rel_ampl[j][j],2./3.);
	    }
	  }
	}
      }
}

double *compute_bispectrum_local(funcparams *pparams, int shape){
  int xi, i, j, k, i1, j1, k1, num, ci, cp;
  int li0, li1, li2, num_x, num_li, np, ni, npd, nid;
  double *red_bispec;

  np = tfvars.num_pol;
  ni = tfvars.num_iso;
  npd = 0; nid = 0;
  for (i = 0; i < np; i++)
    npd += tfvars.pol_arr[shape][i];
  for (i = 0; i < ni; i++)
    nid += tfvars.iso_arr[shape][i];
  li0 = pparams->li[0]; li1 = pparams->li[1]; li2 = pparams->li[2];
  num_x = kintvars.num_x;
  num_li = kintvars.num_li;
  num = npd*npd*npd * nid*nid*(nid+1)/2;
  red_bispec = (double*) malloc(num * sizeof(double));
  for (i = 0; i < num; i++)
    red_bispec[i] = 0.;
  ci = 0;
  for (i1 = 0; i1 < ni; i1++)
    if (tfvars.iso_arr[shape][i1])
      for (j1 = 0; j1 < ni; j1++)
	if (tfvars.iso_arr[shape][j1])
	  for (k1 = j1; k1 < ni; k1++)
	    if (tfvars.iso_arr[shape][k1]){
	      cp = 0;
	      for (i = 0; i < np; i++)
		if (tfvars.pol_arr[shape][i])
		  for (j = 0; j < np; j++)
		    if (tfvars.pol_arr[shape][j])
		      for (k = 0; k < np; k++)
			if (tfvars.pol_arr[shape][k]){
			  for (xi = 0; xi < num_x; xi++){
			    red_bispec[ci*npd*npd*npd + cp] +=
			kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
   		        * (kintvars.beta_iso_arr[j1*np*num_x*num_li + j*num_x*num_li
					     + li1*num_x + xi]
		           * kintvars.beta_iso_arr[k1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]
		           + kintvars.beta_iso_arr[k1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]
		           * kintvars.beta_iso_arr[j1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]);
			    red_bispec[ci*npd*npd*npd + cp] +=
		        kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + j*num_x*num_li
					     + li1*num_x + xi]
		        * (kintvars.beta_iso_arr[j1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
		           * kintvars.beta_iso_arr[k1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]
		           + kintvars.beta_iso_arr[k1*np*num_x*num_li+i*num_x*num_li
					       + li0*num_x + xi]
		           * kintvars.beta_iso_arr[j1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]);
			    red_bispec[ci*npd*npd*npd + cp] +=
		        kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + k*num_x*num_li
					     + li2*num_x + xi]
		        * (kintvars.beta_iso_arr[j1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
		           * kintvars.beta_iso_arr[k1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]
		           + kintvars.beta_iso_arr[k1*np*num_x*num_li+i*num_x*num_li
					       + li0*num_x + xi]
		           * kintvars.beta_iso_arr[j1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]);
			  }
			  cp++;
			}
	      ci++;
	    }
  for (i = 0; i < num; i++)
    red_bispec[i] *= pow(2.0/M_PI, 3);
  // red_bispec[i] *= 2. * pow(2.0/M_PI, 3);
  // This factor of 2 was in older versions of the code taken into
  // account later, but it should be included in red_bispec.
  // Moreover, because of the isocurvature symmetrization, we have to
  // divide by 2 as well, so they cancel.
  return red_bispec;
}

double *compute_bispectrum_local_indep(funcparams *pparams, int shape){
  int xi, i, j, k, i1, j1, k1, num, ci, cp;
  int li0, li1, li2, num_x, num_li, np, ni, npd, nid;
  double *red_bispec;

  np = tfvars.num_pol;
  ni = tfvars.num_iso;
  npd = 0; nid = 0;
  for (i = 0; i < np; i++)
    npd += tfvars.pol_arr[shape][i];
  for (i = 0; i < ni; i++)
    nid += tfvars.iso_arr[shape][i];
  li0 = pparams->li[0]; li1 = pparams->li[1]; li2 = pparams->li[2];
  num_x = kintvars.num_x;
  num_li = kintvars.num_li;
  num = npd*npd*npd * nid*nid*(nid+1)/2;
  red_bispec = (double*) malloc(num * sizeof(double));
  for (i = 0; i < num; i++)
    red_bispec[i] = 0.;
  ci = 0;
  for (i1 = 0; i1 < ni; i1++)
    if (tfvars.iso_arr[shape][i1])
      for (j1 = 0; j1 < ni; j1++)
	if (tfvars.iso_arr[shape][j1])
	  for (k1 = j1; k1 < ni; k1++)
	    if (tfvars.iso_arr[shape][k1]){
	      cp = 0;
	      for (i = 0; i < np; i++)
		if (tfvars.pol_arr[shape][i])
		  for (j = 0; j < np; j++)
		    if (tfvars.pol_arr[shape][j])
		      for (k = 0; k < np; k++)
			if (tfvars.pol_arr[shape][k]){
			  for (xi = 0; xi < num_x; xi++){
			      if (i1==0 && j1==0 && k1==0) {
			    red_bispec[ci*npd*npd*npd + cp] +=
			kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
   		        * (kintvars.beta_adi_arr[j1*np*num_x*num_li + j*num_x*num_li
					     + li1*num_x + xi]
		           * kintvars.beta_adi_arr[k1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]
		           + kintvars.beta_adi_arr[k1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]
		           * kintvars.beta_adi_arr[j1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]);
			    red_bispec[ci*npd*npd*npd + cp] +=
		        kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + j*num_x*num_li
					     + li1*num_x + xi]
		        * (kintvars.beta_adi_arr[j1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
		           * kintvars.beta_adi_arr[k1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]
		           + kintvars.beta_adi_arr[k1*np*num_x*num_li+i*num_x*num_li
					       + li0*num_x + xi]
		           * kintvars.beta_adi_arr[j1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]);
			    red_bispec[ci*npd*npd*npd + cp] +=
		        kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + k*num_x*num_li
					     + li2*num_x + xi]
		        * (kintvars.beta_adi_arr[j1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
		           * kintvars.beta_adi_arr[k1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]
		           + kintvars.beta_adi_arr[k1*np*num_x*num_li+i*num_x*num_li
					       + li0*num_x + xi]
		           * kintvars.beta_adi_arr[j1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]);
                }
                else if (j1==i1 && k1==i1) {
                    red_bispec[ci*npd*npd*npd + cp] +=
			kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
   		        * (kintvars.beta_iso_arr[j1*np*num_x*num_li + j*num_x*num_li
					     + li1*num_x + xi]
		           * kintvars.beta_iso_arr[k1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]
		           + kintvars.beta_iso_arr[k1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]
		           * kintvars.beta_iso_arr[j1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]);
			    red_bispec[ci*npd*npd*npd + cp] +=
		        kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + j*num_x*num_li
					     + li1*num_x + xi]
		        * (kintvars.beta_iso_arr[j1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
		           * kintvars.beta_iso_arr[k1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]
		           + kintvars.beta_iso_arr[k1*np*num_x*num_li+i*num_x*num_li
					       + li0*num_x + xi]
		           * kintvars.beta_iso_arr[j1*np*num_x*num_li+k*num_x*num_li
					       + li2*num_x + xi]);
			    red_bispec[ci*npd*npd*npd + cp] +=
		        kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		        * kintvars.alpha_arr[i1*np*num_x*num_li + k*num_x*num_li
					     + li2*num_x + xi]
		        * (kintvars.beta_iso_arr[j1*np*num_x*num_li + i*num_x*num_li
					     + li0*num_x + xi]
		           * kintvars.beta_iso_arr[k1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]
		           + kintvars.beta_iso_arr[k1*np*num_x*num_li+i*num_x*num_li
					       + li0*num_x + xi]
		           * kintvars.beta_iso_arr[j1*np*num_x*num_li+j*num_x*num_li
					       + li1*num_x + xi]);
                }
                else {red_bispec[ci*npd*npd*npd + cp] = 0;}
			  }
			  cp++;
			}
	      ci++;
	    }
  for (i = 0; i < num; i++)
    red_bispec[i] *= pow(2.0/M_PI, 3);
  // red_bispec[i] *= 2. * pow(2.0/M_PI, 3);
  // This factor of 2 was in older versions of the code taken into
  // account later, but it should be included in red_bispec.
  // Moreover, because of the isocurvature symmetrization, we have to
  // divide by 2 as well, so they cancel.
  return red_bispec;
}

// Adiabatic only, for the moment
double *compute_bispectrum_equil(funcparams *pparams, int shape){
  int xi, i, j, k, num, cp;
  int li0, li1, li2, num_x, num_li, np, npd;
  double *red_bispec;

  np = tfvars.num_pol;
  npd = 0;
  for (i = 0; i < np; i++)
    npd += tfvars.pol_arr[shape][i];
  li0 = pparams->li[0]; li1 = pparams->li[1]; li2 = pparams->li[2];
  num_x = kintvars.num_x;
  num_li = kintvars.num_li;
  num = npd*npd*npd;
  red_bispec = (double*) malloc(num * sizeof(double));
  for (i = 0; i < num; i++)
    red_bispec[i] = 0.;
  cp = 0;
  for (i = 0; i < np; i++)
    if (tfvars.pol_arr[shape][i])
      for (j = 0; j < np; j++)
	if (tfvars.pol_arr[shape][j])
	  for (k = 0; k < np; k++)
	    if (tfvars.pol_arr[shape][k]){
	      for (xi = 0; xi < num_x; xi++){
		red_bispec[cp] +=
		  -1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.alpha_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  -1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.alpha_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  -1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.alpha_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi];

		red_bispec[cp] +=
		  -2. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.delta_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.delta_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.delta_arr[k*num_x*num_li + li2*num_x + xi];

		red_bispec[cp] +=
		  1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.gamma_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.delta_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.gamma_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.delta_arr[j*num_x*num_li + li1*num_x + xi];
		red_bispec[cp] +=
		  1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.gamma_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.delta_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.gamma_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.delta_arr[i*num_x*num_li + li0*num_x + xi];
		red_bispec[cp] +=
		  1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.gamma_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.delta_arr[j*num_x*num_li + li1*num_x + xi];
		red_bispec[cp] +=
		  1. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.gamma_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.delta_arr[i*num_x*num_li + li0*num_x + xi];
	      }
	      cp++;
	    }
  for (i = 0; i < num; i++)
    red_bispec[i] *= 6. * pow(2.0/M_PI, 3);
  return red_bispec;
}

// Adiabatic only, for the moment
double *compute_bispectrum_ortho(funcparams *pparams, int shape){
  int xi, i, j, k, num, cp;
  int li0, li1, li2, num_x, num_li, np, npd;
  double *red_bispec;

  np = tfvars.num_pol;
  npd = 0;
  for (i = 0; i < np; i++)
    npd += tfvars.pol_arr[shape][i];
  li0 = pparams->li[0]; li1 = pparams->li[1]; li2 = pparams->li[2];
  num_x = kintvars.num_x;
  num_li = kintvars.num_li;
  num = npd*npd*npd;
  red_bispec = (double*) malloc(num * sizeof(double));
  for (i = 0; i < num; i++)
    red_bispec[i] = 0.;
  cp = 0;
  for (i = 0; i < np; i++)
    if (tfvars.pol_arr[shape][i])
      for (j = 0; j < np; j++)
	if (tfvars.pol_arr[shape][j])
	  for (k = 0; k < np; k++)
	    if (tfvars.pol_arr[shape][k]){
	      for (xi = 0; xi < num_x; xi++){
		red_bispec[cp] +=
		  -3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.alpha_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  -3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.alpha_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  -3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.alpha_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi];

		red_bispec[cp] +=
		  -8. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.delta_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.delta_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.delta_arr[k*num_x*num_li + li2*num_x + xi];

		red_bispec[cp] +=
		  3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.gamma_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.delta_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.gamma_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.delta_arr[j*num_x*num_li + li1*num_x + xi];
		red_bispec[cp] +=
		  3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.gamma_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.delta_arr[k*num_x*num_li + li2*num_x + xi];
		red_bispec[cp] +=
		  3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.gamma_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.delta_arr[i*num_x*num_li + li0*num_x + xi];
		red_bispec[cp] +=
		  3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.gamma_arr[i*num_x*num_li + li0*num_x + xi]
		  * kintvars.delta_arr[j*num_x*num_li + li1*num_x + xi];
		red_bispec[cp] +=
		  3. * kintvars.w_arr[xi]*kintvars.x_arr[xi]*kintvars.x_arr[xi]
		  * kintvars.beta_iso_arr[k*num_x*num_li + li2*num_x + xi]
		  * kintvars.gamma_arr[j*num_x*num_li + li1*num_x + xi]
		  * kintvars.delta_arr[i*num_x*num_li + li0*num_x + xi];
	      }
	      cp++;
	    }
  for (i = 0; i < num; i++)
    red_bispec[i] *= 6. * pow(2.0/M_PI, 3);
  return red_bispec;
}

double *compute_cl(funcparams *pparams){
  int count, h, i, j, k;
  int num_pol = tfvars.num_pol;
  int num_iso = tfvars.num_iso;
  int np = tfvars.num_pol_do;
  int ni = tfvars.num_iso_do;
  double *cl, amplitude;

  cl = (double*) malloc(np*(np+1)/2 * ni*(ni+1)/2 * sizeof(double));
  count = 0;
  for (h = 0; h < num_iso; h++)
    if (tfvars.iso_summ[h]){
      pparams->iso[0] = h;
      for (i = h; i < num_iso; i++)
	if (tfvars.iso_summ[i]){
	  pparams->iso[1] = i;
	  pparams->power = -1. + tfvars.red_tilt[h][i];
	  if (h==0 && i==0) {
        amplitude = 1.;
      } else if (h==0 || i==0) {
        amplitude = tfvars.rel_ampl_ai;
      } else {
        amplitude = tfvars.rel_ampl_iso;
	  }
	  for (j = 0; j < num_pol; j++)
	    if (tfvars.pol_summ[j]){
	      pparams->pol[0] = j;
	      for (k = j; k < num_pol; k++)
		if (tfvars.pol_summ[k]){
		  pparams->pol[1] = k;
		  cl[count] = (2.0/M_PI) * tfvars.rel_ampl[h][i] * pow(tfvars.pivot,-1) * amplitude
		              * int_1d(&tf_integrand, pparams->qmin,
				       pparams->qmax, PTS_CL, NSUB_CL, pparams);
		  count++;
		}
	    }
	}
    }
  return cl;
}
